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We study the transition from a Mott insulator to a superfluid in both the two- and the three- 
dimensional Bose-Hubbard model at zero temperature, employing the method of the effective po- 
tential. Converting Kato's perturbation series into an algorithm capable of reaching high orders, we 
obtain accurate critical parameters for any integer filling factor. Our technique allows us to monitor 
both the approach to the mean-field limit by considering spatial dimensionalities d > 3, and to the 
quantum rotor limit of high filling, which refers to an array of Josephson junctions. 
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The Bose-Hubbard model, describing interacting Bose 
particles moving on a tight-binding lattice, has drawn 
much attention, especially after its experimental realiza- 
tion with ultracold bosonic atoms in optical potentials 
(see Ref. [ij and references therein) . This clean defectless 
setup, which allows for precise control of its parameters, 
has opened up new testing ground for quantum many- 
body physics. The pure Bose-Hubbard system reflects 
the competition between the potential energy due to the 
repulsive on-site interaction among the Bosons, which 
tends to suppress density fluctuations and to localize the 
particles, and the kinetic energy associated with tunnel- 
ing processes between neighboring lattice sites, which try 
to delocalize the particles and to reduce phase fluctua- 
tions. Denoting the on-site interaction energy of a pair 
of particles sitting at the same site by U, and the hop- 
ping matrix element by J, the model's grand canonical 
Hamiltonian is written in dimcnsionlcss form as 



Hbh — 



i J2 n,{n, - 1) - ti/UY^h, -J/U J2 
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(1) 

where indices label the sites of a d-dimensional lattice, 
which we take as hypercubic, and the sum over (j, k) ex- 
tends over nearest neighbors. As usual, a| and a^- are the 
creation and annihilation operators for a Boson at site j, 
and fij = OjUj is the number operator at that site. The 
chemical potential fi here is site- independent. At zero 
temperature one finds a series of Mott phases at suffi- 
ciently small values of J/U, characterized by a fixed fill- 
ing of an integer number of particles per site, depending 
on the value of n/U. A Mott state has zero compress- 
ibility, due to an energy gap separating the ground state 
from the particle and hole excitations, so that it costs 
energy to move a particle through the system. Upon in- 
crasing the ratio J/U , the competition between potential 
and kinetic energy leads to a quantum phase transition: 



At the phase boundary (J/[/)pb the gap closes, so that 
particle delocalization becomes favorable, and the system 
Bose-condenses into a superfluid state for d > 2 In 
optical lattice experiments performed so far, this transi- 
tion has been induced by varying the lattice depth Q , as 
in the pioneering work by Greiner et al. 0] , and by shak- 
ing the lattice periodically in time with slowly varying 
amplitude Q, as done recently by Zenesini et al. 0] 

Despite the apparent simplicity of the Hamiltonian ([1]) , 
a precise calculation of its critical parameters for different 
dimensionalities d and filling factors g poses severe chal- 
lenges, so that the determination of the phase diagram 
in the J/U-^/U-plane has become a major benchmark 
problem for computational many-body physics. Recent 
quantum Monte Carlo (QMC) simulations have yielded 
critical parameters with record accuracy for 5 = 1 0i Si ■ 
A previous strong-coupling expansion had led to reliable 
analytical results to third order in J/U 9], and later 
was extended to higher orders in one and two dimensions 
for g = I and g = 2 [iq|. Techniques using the density- 
matrix renormalization group (DMRG) allow one to treat 
fairly large systems in one dimension jlll . 12 , 1^ , but up 
to now have remained restricted to low filling. So far, ac- 
curate critical data for the three-dimensional (3D) system 
with experimentally relevant higher filling factors g > 1 
have remained particularly hard to obtain. 

In this contribution we show that a specific adaption 
of high-order many-body perturbation theor y, b ased on 



Kato's formulation of the perturbation series [14|, |15l| and 
using the concept of the order parameter, enables one to 
investigate Bose-Hubbard systems with arbitrary integer 
filling factor. In principle, the technique is applicable to 
any type of lattice, in any dimension. We first briefly 
sketch the method, and present our results for both 2D 
and 3D lattices. We then numerically monitor the ap- 
proach to the mean-field limit of high lattice dimension, 
and to the quantum rotor limit of high filling 
which describes a Josephson junction array [18|. 
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Our starting point is the method of the effective po- 
tential 191, as considered recently by dos Santos and 
Pelster [20|. Adding source terms to the Bose-Hubbard 
Hamiltonian ([1]) which attempt to add particles with uni- 
form strength x to each site, or to remove them with 
strength %* according to 



H^YiiXi X*) = Ho + Htun + (x*aj + xa 



(2) 



then expanding the grand-canonical free energy F — 
{Hbk) at zero temperature into a power series in the 
hopping parameter J/U and the sources Xi X* i one has 



FiJ/U, x,X*) = M FoiJ/U) + J2 C2niJ/U)\x 

\ n 

for a lattice of M sites, with coefficients 



\2n 



(3) 



(4) 



The order parameter ip now specifies the change of F in 
response to a variation of the sources, 
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(5) 



M dx* 

while the effective potential T — F/M — ilj*x — V'X* is the 
Legendre transform of F, with ip and ip* as independent 
variables. With the help of Eqs. (O and ([3]) one gets the 
familiar Landau form 

1 



C2 



C4 



(6) 



Since dT/dtp — — x* and dT/dip* ~ — and since the 
original Bose-Hubbard system is recovered by setting 
X = X* = 0, the system adopts that order parameter 
which minimizes F. Unless fj,/U is integer, one has C2 < 
for sufficiently small J/U, whereas C4 > 0, so that one 
finds a Mott regime with ■0 = 0. Upon increasing J/U, 
the system enters the superfluid phase when acquires 
a nonzero value, indicating long-range phase coherence. 
Hence, the phase boundary is determined by that J/U 
for which the minimum of the expression ([6]) starts to 
deviate from \tp\'^ — 0, which occurs when the coefficient 
— l/c2 of vanishes. In effect, one has to identify that 
scaled hopping strength J/U for which the susceptibility 
C2 divergesjthis divergence marks the quantum phase 
transition [20|. 

For computing C2 we resort to Kato's formulation of 
the perturbation series 14, [l^, starting from the site- 
diagonal Hamiltonian Hq. For integer filling factor g, 
its ground state |m) is a product of local Fock states 
with g particles sitting at each site. In general, when the 
system is subjected to some perturbation V, the nth- 



order correction to its energy is given by the trace [14 1 



Sf"^, = tr 

|m) 



(7) 



where the sum runs over all possible sets of nonnegative 
integers ai which obey ae = n—1. The operators S" 
are defined by 



E 



-|m)(m| 

{En, - E,Y 



for a = 
for a > 



(8) 



with and Ei denoting the unperturbed energies of 
the ifo-eigenstates |m) and respectively. This ex- 
pression ([7]) can be understood as a sum over chains of 
processes mediated by the operators V . Each process 
chain leads from the Mott-insulator state |m) over dif- 
ferent intermediate states \i) back to |m). Such chains 
can be represented by abstract diagrams, with only con- 
nected diagrams contributing to the sum, as stated by the 
linked-cluster theorem 21|. Each diagram has a certain 
weight depending on the lattice's type and dimensional- 
ity. For example, diagrams for the energy correction due 
to tunneling consist merely of closed loops of individual 
tunneling processes. In contrast, for calculating C2 the 
augmented Hamiltonian ([2]) prompts us to set 



E 



X + xa) 



(9) 



V = -J/U J2 «]' 

{3,k) j 

Because we are aiming at the coefficient of |xP in Eq. 
we only need to take into account terms containing ex- 
actly one creation and one annihilation process. This se- 
lection yields C2{J/U) = J^i, ^2^\j /UY as a series in the 
tunneling parameter J /U . The only relevant third-order 
diagram thus consists of one creation of a Boson (•), one 
tunneling process (— >), and one annihilation (x). The 
fourth-order diagrams then read 

• X , •x , (10) 

with the second diagram indicating chains for which cre- 
ation and annihilation take place at the same site. The 
computational effort increases quickly with the order: 
For = 8, say, all permutations of up to ten different 
processes (8 — >, 1 •, 1 x) encoded in the diagrams have 
to be evaluated. 

An instructive example for illustrating this scheme oc- 
curs in the limit of infinite lattice dimensionality d. Here 
the diagrams containing "back and forth" tunneling pro- 
cesses (analogous to the second diagram pO)) ) do not con- 
tribute to the sum, because they acquire vanishing weight 
for d 00. The remaining diagrams simply are 

• X , • — > X , • X , . . . , • . . . ^ X . (11) 

Being one-particle reducible, they factorize into their 
one-particle irreducible contributions (2^, [l^l : 

\2 



X = (-1) (.x)^ 
X = i-lYi'xf 



(12) 



.(^)'^x = (-l)''(.x) 
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FIG. 1: Logarithm of the coefficients — aj"' for filling factors 
g — 1, 10, 50 in two and three dimensions, with linear fits. 
The chemical potential is chosen as /i/U = g — 0.5. 



For each tunneling process one has an additional fac- 
tor 2d, since there exist 2c? directions on a d-dimensional 
rectangular lattice. The resulting series for C2{J /U) is ge- 



ometric, because a. 



{u-l) 



l/(2da2°^) is constant; 



this ratio determines its radius of convergence and hence 
directly gives the phase boundary: 



2d{J/U\,^ = 



(g-A^/t/)W£/-ff+l) 

^Ji/u + l 



(13) 



which is precisely the mean- field result 0, [l3| ■ 

We have devised an algorithm for efficiently generating 
and evaluating all diagrams up to some order for any lat- 
tice dimension d. In two and three dimensions we obtain 
(negative) coefficients 0^2^ which form almost perfect ge- 
ometric series, as depicted in Fig.[T]for 5 = 1, 10, and 50. 



If the ratio a. 



(y-l) 



were constant, it would equal the 



phase boundary as in the example above. But since now 
this ratio changes slightly with the number v of tunneling 
processes taken into account, we carry out an extrapola- 
tion over l/i/ by making a linear fit based on the orders 
1 to 8 in J /U (3 to 10 in V), as illustrated by the central 
inset in Fig.[2l Different selections of the orders employed 
(e.g., 2 to 8 in J /U) lead to very similar results, with an 
uncertainty of about 1% in 3D, and 2% in 2D. The main 
part of Fig. [2] shows the phase boundary thus obtained 
for the 3D case at unit filling, together with some approx- 
imants for finite orders. The tip of the lobe corresponds 
to the critical parameter {J/U)c, for which QMC calcu- 
lations have provided a highly accurate reference value: 
{ J/U)c = 0.03408(2) for 5 = 1 0. Our data match this 
value fairly well, as emphasized by the lower right inset. 

Critical parameters obtained for higher filling g in two 
and three dimensions are collected in Tab. HI With in- 
creasing g, the critical chemical potential {^/U)c ap- 
proaches g — 0.5, due to the fact that there is exact 
particle- hole symmetry for 5 — > 00. Some correspond- 
ing Mott lobes are depicted in Fig. [3l for 5 = 1, QMC 
data 0, S] are included for comparison. 

Our technique permits us to reach higher dimension- 
alities d > 3, thus uncovering how the mean-field limit 
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FIG. 2: Phase boundary for the 3D model with unit filling, 
as determined from the ratios q'''~^-'/q'''' for finite orders u, 
together with the extrapolation to = 00 (extr). The in- 
set at the right bottom magnifies the tip of the lobe, and 
demonstrates the convergence to the QMC result 0] (dashed 
vertical line). The central inset illustrates the extrapolation 
of a'^-^V"^"' to (J/(7)c for d = 2 (upper data) and d = 3 
(lower data). Observe that the data for d = 3 fluctuate less. 



TABLE L Critical values {fi/U)c and {J/U)c for various filling 
factors g. For locating the tip of the respective Mott lobe, n/U 
has been varied in steps of 0.001. Relative errors of {J/U)c 
are less than 1% for d = 3, and less than 2% for d — 2. 





d = 2 


d = 3 
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(J/U). 




iJ/u). 


1 


0.376 


5.909E-002 


0.393 


3.407E-002 


2 


1.427 


3.480E-002 


1.437 


2.007E-002 


3 


2.448 


2.473E-002 


2.455 


1.427E-002 


4 


3.460 


1.920E-002 


3.465 


1.108E-002 


5 


4.470 


1.569E-002 


4.472 


9.055E-003 


10 


9.483 


8.208E-003 


9.485 


4.736E-003 


20 


19.491 


4.202E-003 


19.492 


2.425E-003 


50 


49.496 


1.706E-003 


49.497 


9.842E-004 


100 


99.498 


8.571E-004 


99.498 


4.946E-004 


1000 


999.50 


8.609E-005 


999.50 


4.968E-005 


10000 


9999.50 


8.613E-006 


9999.50 


4.970E-006 



is approached, and high filling factors g 3> 1. In the 
latter regime, the phases at the individual sites become 
well defined, so that the Bose-Hubbard model reduces 
to a quantum rotor model containing a single parameter 
gJ/U, and describing a Josephson junction array [lR[l^. 
Figure m indeed reveals that the products 2dg{J/U)c re- 
main almost constant when g exceeds 100, with limiting 
values 0.345 for d = 2 and 0.299 for d — 3 falling sig- 
nificantly above the mean-field prediction of 1/4, which 
follows from Eq. (I13p . Even for d — 10, the data still 
exceed the mean- field result by 4%. 

To conclude, diagrammatic many-body perturbation 
theory based on Kato's series ([7]), though impractical to 
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gJ/U 

FIG. 3: Mott lobes for d — 2 (upper panel) and d — 3 (below) 
with various g. Dashed lines mark the quantum rotor limit 
{fJ,/U)c ~ g — 0.5 of the critical chemical potential. The lobes' 
tips are magnified in the inset, illustratiriE the convergence of 
g{J/U)c. For unit filling, QMC data 0,|1| are included. 
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FIG. 4: Critical product 2dg{J/U)c for d = 2, 3, 5, and 10 
vs. g, together with the mean-field limit. Even for d = 10, the 
large-y-Umit still exceeds the mean-field prediction by 4%. 



work out analytically in high orders, becomes a powerful 
and accurate tool when turned into a numerically exe- 
cutable algorithm. The merit of this technique rests not 
only in the fact that it enables one to access regimes 
which could not be reached before, such as experimen- 
tally important filling factors g > 1 or the crossover 
to the quantum rotor dynamics depicted in Fig. |4l but 
also in its great flexibility. For instance, with appro- 



priately constructed diagrams it also yields correlation 
functions. Thus, the applicability of this approach is by 
no means exhausted by the present calculation of the 
Bose-Hubbard phase diagram. 
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